Debugging

Import Data

corrcounts_merge <- readRDS("~/VersionControl/senescence_benchmarking/Data/corrcounts_merge.rds")
metadata_merge <- readRDS("~/VersionControl/senescence_benchmarking/Data/metadata_merge.rds")
SenescenceSignatures <- readRDS("~/VersionControl/senescence_benchmarking/CommonFiles/SenescenceSignatures_divided_newCellAge.RDS")

library(markeR)
library(ggplot2)
library(ggpubr)
library(edgeR)

Loading required package: limma

?markeR

ℹ Rendering development documentation for "markeR"

Scores

?CalculateScores

ℹ Rendering development documentation for "CalculateScores"

df_logmedian <- CalculateScores(data = corrcounts_merge, metadata = metadata_merge, method = "logmedian", gene_sets = SenescenceSignatures)

senescence_triggers_colors <- c(
  "none" = "#E57373",  # Soft red  
  "Radiation" = "#BDBDBD",  # Medium gray  
  "DNA damage" = "#64B5F6",  # Brighter blue  
  "Telomere shortening" = "#4FC3F7",  # Vivid sky blue  
  "DNA demethylation" = "#BA68C8",  # Rich lavender  
  "Oxidative stress" = "#FDD835",  # Strong yellow  
  "Conditioned Medium" = "#F2994A",  # Warm orange  
  "Oncogene" = "#81C784",  # Medium green  
  "Lipid Accumulation" = "#E57373",  # Coral  
  "Calcium influx" = "#26A69A",  # Deep teal  
  "Plasma membrane dysruption" = "#D32F2F",  # Strong salmon  
  "OSKM factors" = "#FFB74D",  # Bright peach  
  "YAP KO" = "#9575CD"  # Deep pastel purple  
)

cellTypes_colors <- c(
  "Fibroblast" = "#FF6961",   # Strong Pastel Red  
  "Keratinocyte" = "#FFB347", # Strong Pastel Orange  
  "Melanocyte" = "#FFD700",   # Strong Pastel Yellow  
  "Endothelial" = "#77DD77",  # Strong Pastel Green  
  "Neuronal" = "#779ECB",     # Strong Pastel Blue  
  "Mesenchymal" = "#C27BA0"   # Strong Pastel Purple  
)

cond_cohend <- list(A=c("Senescent"), # if no variable is defined, will be the first that appears in the ggplot
                    B=c("Proliferative","Quiescent"))

PlotScores(ResultsList = df_logmedian, ColorVariable = "CellType", GroupingVariable="Condition",  method ="logmedian", ColorValues = cellTypes_colors, ConnectGroups=TRUE, ncol = 6, nrow = 2, widthTitle=20, y_limits = NULL, legend_nrow = 2,xlab=NULL, cond_cohend = cond_cohend)

Try scores with bidirectional signatures

bidirectsigs <- readRDS("~/VersionControl/senescence_benchmarking/CommonFiles/SenescenceSignatures_complete_newCellAge.RDS")
for (sig in names(bidirectsigs)){
  sigdf <- bidirectsigs[[sig]]
  sigdf <- sigdf[,1:2] # remove the third column, if applicable
  if(any(sigdf[,2]=="not_reported")){
    sigdf <- sigdf[,1]
    bidirectsigs[[sig]] <- sigdf
    next 
  }
  sigdf[,2] <- ifelse(sigdf[,2]=="enriched",1,-1)
  bidirectsigs[[sig]] <- sigdf
}
bidirectsigs

$CellAge

$CSgene
  [1] "TP53"       "TERF2"      "MAPK14"     "CDKN2A"     "CDKN1A"     "CCNE1"      "CCNA1"      "MAPKAPK5"   "CBX4"       "TXN"       
 [11] "TBX2"       "STAT3"      "SRF"        "BMI1"       "MAP2K4"     "MAP2K6"     "MAP2K3"     "MAPK8"      "MAPK3"      "MAPK1"     
 [21] "PRKCD"      "PML"        "OPA1"       "ATM"        "MDM2"       "CXCL8"      "IL6"        "IGFBP7"     "ID1"        "HRAS"      
 [31] "H2AFX"      "POT1"       "SIRT1"      "KDM6B"      "PLA2R1"     "EZH2"       "E2F3"       "E2F1"       "CEBPB"      "CDKN2D"    
 [41] "CDKN2B"     "CDKN1B"     "CDK6"       "CDK4"       "CDK2"       "CDC42"      "RBX1"       "CDC27"      "CDK1"       "MAML1"     
 [51] "CD44"       "MAD2L1BP"   "MAP4K4"     "AIM2"       "RECQL4"     "ARHGAP18"   "KL"         "MAPKAPK2"   "AURKB"      "SLC16A7"   
 [61] "CCNE2"      "HIST1H2BJ"  "HIST1H3F"   "CCNA2"      "MCM3AP"     "CDC16"      "TSC22D1"    "CBS"        "TNFSF13"    "CTNNAL1"   
 [71] "EED"        "PNPT1"      "CDC23"      "RNASET2"    "TP63"       "CAV1"       "MKNK1"      "TSLP"       "HIST1H2BK"  "PPM1D"     
 [81] "HAVCR2"     "CBX2"       "KDM2B"      "DPY30"      "C2orf40"    "YPEL3"      "HIST2H4A"   "HIST1H4L"   "HIST1H4E"   "HIST1H4B"  
 [91] "HIST1H4H"   "HIST1H4C"   "HIST1H4J"   "HIST1H4K"   "HIST1H4F"   "HIST1H4D"   "HIST1H4A"   "HIST1H3B"   "HIST1H3H"   "HIST1H3J"  
[101] "HIST1H3G"   "HIST1H3I"   "HIST1H3E"   "HIST1H3C"   "HIST1H3D"   "HIST1H3A"   "HIST2H2BE"  "HIST1H2BO"  "HIST1H2BC"  "HIST1H2BI" 
[111] "HIST1H2BH"  "HIST1H2BE"  "HIST1H2BF"  "HIST1H2BM"  "HIST1H2BN"  "HIST1H2BL"  "HIST1H2BG"  "HIST2H2AC"  "HIST2H2AA3" "HIST1H2AB" 
[121] "HIST1H2AC"  "HIST1H2AJ"  "HIST1H4I"   "HIST3H3"    "CALR"       "HMGA2"      "PHC3"       "KAT6A"      "EHMT1"      "SMC6"      
[131] "AIMP2"      "CALCA"      "DEK"        "MAPKAPK3"   "ZNF148"     "YY1"        "WRN"        "WNT5A"      "NR1H2"      "UBE3A"     
[141] "UBE2E1"     "UBE2D1"     "UBC"        "UBB"        "UBA52"      "CDC26P1"    "TYMS"       "TWIST1"     "HIRA"       "RPS27AP11" 
[151] "HIST2H2AA4" "TP73"       "TOP 1,00"   "TNF"        "TGFB2"      "TGFB1"      "TFDP1"      "TERT"       "TERF1"      "BUB1B"     
[161] "BUB1"       "TCF3"       "TBX3"       "TAGLN"      "STAT6"      "STAT1"      "BRAF"       "SREBF1"     "BRCA1"      "SP1"       
[171] "SOX5"       "SOD2"       "SNAI1"      "SMARCB1"    "SMARCA2"    "HIST2H3D"   "PHC1P1"     "ACD"        "SKIL"       "LOC649620" 
[181] "SLC13A3"    "LOC647654"  "SMURF2"     "ANAPC1"     "SHC1"       "CPEB1"      "H3F3AP6"    "ZMAT3"      "RBBP4P1"    "SRSF3"     
[191] "SRSF1"      "SATB1"      "S100A6"     "RXRB"       "RRM2"       "RRM1"       "RPS27A"     "RPS6KA3"    "RPS6KA2"    "RPS6KA1"   
[201] "RPL5"       "RNF2"       "RIT1"       "RING1"      "BCL2L1"     "RELA"       "BCL2"       "CCND1"      "RBP2"       "RBL2"      
[211] "RBL1"       "RBBP7"      "RBBP4"      "NTN4"       "RB1"        "IL21"       "RAN"        "RAF1"       wrap_title_aux "TNRC6C"    
[221] "KIAA1524"   "EP400"      "CNOT6"      "CBX8"       "PTEN"       "SEPN1"      "BACH1"      "PSMB5"      "PROX1"      "PRL"       
[231] "MAP2K7"     "MAP2K1"     "MAPK10"     "MAPK9"      "MAPK11"     "MAPK7"      "PRKDC"      "RNF114"     "PRKCI"      "ATF7IP"    
[241] "MFN1"       "PRKAA2"     "CDKN2AIP"   "RBM38"      "PRG2"       "HIST2H4B"   "HJURP"      "TMEM140"    "PBRM1"      "Mar-05"    
[251] "PPARG"      "PPARD"      "POU2F1"     "TERF2IP"    "ERRFI1"     "H2BFS"      "PLK1"       "PLAUR"      "PIN1"       "PIM1"      
[261] "PIK3CA"     "PHB"        "PGR"        "PGD"        "PIAS4"      "PDGFB"      "SIRT6"      "ANAPC11"    "ANAPC7"     "ANAPC5"    
[271] "WNT16"      "FZR1"       "ZBTB7A"     "ERGIC2"     "PCNA"       "FIS1"       "PAX3"       "NOX4"       "MINK1"      "PEBP1"     
[281] "YBX1"       "NINJ1"      "NFKB1"      "H2AFB1"     "NDN"        "NCAM1"      "NBN"        "MYC"        "MYBL2"      "MSN"       
[291] "ASS1"       "LOC441488"  "MRE11A"     "MOV10"      "MMP7"       "MIF"        "MAP3K5"     "MAP3K1"     "MECP2"      "MCL1"      
[301] "MAGEA2"     "SMAD9"      "SMAD7"      "SMAD6"      "SMAD5"      "SMAD4"      "SMAD3"      "SMAD2"      "SMAD1"      "MAD2L1"    
[311] "MXD1"       "MIR34A"     "MIR30A"     "MIR299"     "MIR29A"     "MIR22"      "MIR217"     "MIR21"      "MIR205"     "MIR203A"   
[321] "MIR191"     "MIR146A"    "MIR141"     "MIR10B"     "ARNTL"      "LMNB1"      "LMNA"       "LGALS9"     "RHOA"       "KRT5"      
[331] "KRAS"       "KIT"        "KIR2DL4"    "KCNJ12"     "JUN"        "JAK2"       "ITGB4"      "IRS1"       "IRF7"       "IRF5"      
[341] "IRF3"       "ING1"       "IDO1"       "ILF3"       "IL15"       "IL12B"      "CXCR2"      "IL4"        "IGFBP5"     "IGFBP3"    
[351] "IGFBP1"     "IGF1R"      "IGF1"       "H3F3AP5"    "IFNG"       "IFI16"      "IDH1"       "ID2"        "HIST2H3A"   "BIRC5"     
[361] "HSPB1"      "HSPA9"      "HSPA5"      "HSPA1A"     "APEX1"      "HNRNPA1"    "FOXA3"      "FOXA2"      "FOXA1"      "HMGA1"     
[371] "HIF1A"      "ANXA5"      "HELLS"      "HDAC1"      "H3F3B"      "H3F3A"      "HIST1H2BB"  "HIST1H2BD"  "H2AFZ"      "HIST1H2AD" 
[381] "HIST1H2AE"  "ANAPC4"     "ANAPC2"     "UBN1"       "SENP1"      "GUCY2C"     "GSK3B"      "UHRF1"      "BRD7"       "NSMCE2"    
[391] "PTRF"       "GPI"        "GNAO1"      "RPS6KA6"    "TNRC6A"     "AGO2"       "B3GAT1"     "DNAJC2"     "GJA1"       "AGO1"      
[401] "EHF"        "TINF2"      "LDLRAP1"    "ULK3"       "GAPDH"      "ABI3BP"     "ASF1A"      "HIST1H2BA"  "G6PD"       "ACKR1"     
[411] "MTOR"       "CDC26"      "CNOT6L"     "FOS"        "CABIN1"     "MORC3"      "SUZ12"      "NPTXR"      "CBX6"       "SIRT3"     
[421] "CRTC1"      "PPP1R13B"   "SUN1"       "SMC5"       "TNRC6B"     "FOXO1"      "FOXM1"      "TNIK"       "SCMH1"      "DKK 1,00"  
[431] "FGFR2"      "FGF2"       "HEPACAM"    "FANCD2"     "EWSR1"      "ETS2"       "ETS1"       "ESR2"       "ERF"        "AKT1"      
[441] "EREG"       "ERBB2"      "ENG"        "ELN"        "CRTC2"      "EIF5A"      "EGR1"       "EGFR"       "EEF1B2"     "AGO4"      
[451] "AGO3"       "EEF1A1"     "PHC2"       "PHC1"       "ABCA1"      "E2F2"       "DUSP6"      "DUSP4"      "HBEGF"      "AGT"       
[461] "DNMT3A"     "AGER"       "DKC1"       "DAXX"       "CYP3A4"     "CTSZ"       "CTSD"       "CSNK2A1"    "E2F7"       "PARP1"     
[471] "HIST3H2BB"  "HIST2H3C"   "JDP2"       "HIST4H4"    "CLU"        "CKB"        "RASSF1"     "CHEK1"      "TOPBP1"     "UBE2C"     
[481] "KIF2C"      "BTG3"       "EHMT2"      "GADD45G"    "NEK6"       "ZMYND11"    "SPINT2"     "CENPA"      "AGR2"       "CEBPG"     
[491] "HYOU1"      "TADA3"      "MCRS1"      "NDRG1"      "ANAPC10"    "CDKN2C"     "ZMPSTE24"   "PSMD14"     "NAMPT"      "RAD50"     
[501] "TRIM10"     "DNM1L"      "BCL2L11"   

$GOBP_CELLULAR_SENESCENCE
  [1] "AKT3"     "MIR543"   "CDK2"     "CDK6"     "CDKN1A"   "ZMPSTE24" "CDKN1B"   "CDKN2A"   "CDKN2B"   "CITED2"   "KAT5"     "PLK2"    
 [13] "NEK6"     "ZNF277"   "CGAS"     "COMP"     "MAPK14"   "VASH1"    "PLA2R1"   "SMC5"     "SIRT1"    "MORC3"    "NUP62"    "ABL1"    
 [25] "ULK3"     "RSL1D1"   "FBXO5"    "FBXO4"    "MAGEA2B"  "NSMCE2"   "H2AX"     "HLA-G"    "HMGA1"    "HRAS"     "ID2"      "IGF1R"   
 [37] "ING2"     "KIR2DL4"  "ARG2"     "LMNA"     "BMAL1"    "MIR10A"   "MIR146A"  "MIR17"    "MIR188"   "MIR217"   "MIR22"    "MIR34A"  
 [49] "MAGEA2"   "MAP3K3"   "MAP3K5"   "MIF"      "MNT"      "ATM"      "NPM1"     "YBX1"     "OPA1"     "PAWR"     "ABI3"     "FZR1"    
 [61] "WNT16"    "SIRT6"    "PML"      "PRMT6"    "PRELP"    "PRKCD"    "MAPK8"    "MAPK11"   "MAPK9"    "MAPK10"   "MAP2K1"   "MAP2K3"  
 [73] "MAP2K6"   "MAP2K7"   "B2M"      "ZMIZ1"    "PTEN"     "MIR20B"   "RBL1"     "BCL6"     "MAP2K4"   "BMPR1A"   "SPI1"     "SRF"     
 [85] "BRCA2"    "NEK4"     "TBX2"     "TBX3"     "MIR590"   "TERC"     "TERF2"    "TERT"     "TOP2B"    "TP53"     "TWIST1"   "WNT1"    
 [97] "WRN"      "SMC6"     "KAT6A"    "ZKSCAN3"  "HMGA2"    "CALR"     "YPEL3"    "ECRG4"    "MAPKAPK5" "TP63"     "PNPT1"    "DNAJA3"  
[109] "EEF1E1"   "NUAK1"   

$GOBP_NEGATIVE_REGULATION_OF_CELLULAR_SENESCENCE

$GOBP_POSITIVE_REGULATION_OF_CELLULAR_SENESCENCE

$HernandezSegura

$REACTOME_CELLULAR_SENESCENCE
  [1] "CDC27"    "E2F2"     "SCMH1"    "MRE11"    "MAP2K3"   "MAPK9"    "ANAPC4"   "MAP2K4"   "MAP4K4"   "RPS6KA2"  "UBE2D1"   "EED"     
 [13] "MAP2K7"   "TNRC6C"   "MAPKAPK5" "ANAPC5"   "TNRC6A"   "TINF2"    "AGO1"     "CDC23"    "CABIN1"   "MAPK1"    "HIRA"     "TNRC6B"  
 [25] "E2F1"     "RBBP7"    "MAPK3"    "ACD"      "NBN"      "CCNE1"    "FZR1"     "ERF"      "CDK6"     "H2AZ2"    "EZH2"     "MAPK8"   
 [37] "UBE2S"    "MAP2K6"   "NFKB1"    "MAPK10"   "ANAPC15"  "CDKN1B"   "PHC1"     "ASF1A"    "MAPK14"   "E2F3"     "LMNB1"    "RAD50"   
 [49] "TFDP2"    "MAPKAPK3" "IL1A"     "RPS6KA1"  "UBN1"     "RNF2"     "CDKN2C"   "CDK2"     "H1-3"     "H1-1"     "H2BC11"   "CDKN1A"  
 [61] "ID1"      "AGO3"     "POT1"     "CDKN2D"   "CDC16"    "H3-3B"    "KDM6B"    "TERF2"    "CCNA1"    "PHC2"     "AGO4"     "ETS1"    
 [73] "CDK4"     "MDM2"     "IL6"      "TXN"      "HMGA1"    "RB1"      "MINK1"    "TP53"     "ANAPC11"  "CBX8"     "CBX4"     "RPS27A"  
 [85] "CCNA2"    "H2BC1"    "TERF1"    "CDKN2B"   "CDKN2A"   "ATM"      "HMGA2"    "UBC"      "VENTX"    "ANAPC1"   "TNIK"     "MOV10"   
 [97] "ETS2"     "H2BC5"    "H4C8"     "RBBP4"    "MAPKAPK2" "H3-3A"    "IGFBP7"   "ANAPC10"  "ANAPC16"  "MAPK7"    "TERF2IP"  "H3-4"    
[109] "BMI1"     "H1-4"     "STAT3"    "CXCL8"    "UBE2E1"   "UBB"      "FOS"      "IFNB1"    "CEBPB"    "KAT5"     "RELA"     "PHC3"    
[121] "CBX2"     "UBE2C"    "CCNE2"    "ANAPC2"   "CDC26"    "RPS6KA3"  "JUN"      "SUZ12"    "H2AC6"    "H2BC4"    "EHMT1"    "EP400"   
[133] "H3C13"    "CBX6"     "H2AC20"   "H1-5"     "H2BC21"   "H2BC13"   "MAPK11"   "SP1"      "H1-2"     "H2AX"     "H1-0"     "ANAPC7"  
[145] "H2AC7"    "H2BC26"   "H4C3"     "H3C12"    "H4C11"    "H3C4"     "MAP3K5"   "H4C16"    "H2BC12"   "TFDP1"    "MDM4"     "H3C14"   
[157] "H3C15"    "RING1"    "EHMT2"    "UBA52"    "H2AJ"     "H4C15"    "H4C14"    "H4C12"    "H2BC14"   "H2BC8"    "H3C8"     "H2AB1"   
[169] "H2BC6"    "H4C6"     "H2BC17"   "H3C6"     "H4C13"    "H3C11"    "H2BC9"    "H3C1"     "H4C9"     "H2AC14"   "H2BC3"    "H4C5"    
[181] "H2AC8"    "H4C4"     "H2BC7"    "H3C7"     "H2AC4"    "H2BC10"   "H4C1"     "H4C2"     "H3C10"    "MIR24-2"  "MIR24-1"  "H3C2"    
[193] "H3C3"     "H2AC18"   "H2AC19"  

$SAUL_SEN_MAYO

$SeneQuest

NA

wrap_title_aux <- function(title, width = 30) {
  if (nchar(title) <= width) {
    return(title)  # No need to wrap if it fits
  }

  wrapped_title <- ""
  while (nchar(title) > width) {
    # Find positions of capital letters and symbols near the wrap point
    capital_pos <- gregexpr("[A-Z]", title)[[1]]
    symbol_pos <- gregexpr("(_|-|:)", title)[[1]]

    # Check for symbol breaks within the last few characters (width - 5 to width)
    valid_symbol_breaks <- symbol_pos[symbol_pos >= (width - 5) & symbol_pos <= width]

    if (length(valid_symbol_breaks) > 0) {
      # If a suitable symbol is found, break at the first valid symbol
      break_at <- valid_symbol_breaks[1]
    } else {
      # If no suitable symbol, look for capital letters within the same range
      valid_capital_breaks <- capital_pos[capital_pos >= (width - 5) & capital_pos <= width]

      if (length(valid_capital_breaks) > 0) {
        # If a capital letter is found, break just before the capital letter
        break_at <- valid_capital_breaks[1] - 1
      } else {
        # If no suitable symbol or capital letter, break at width
        break_at <- width
      }
    }

    # Append the wrapped line
    wrapped_title <- paste0(wrapped_title, substr(title, 1, break_at), "\n")

    # Update title with the remaining text after the break
    title <- substr(title, break_at + 1, nchar(title))
  }

  # Add the remaining part of the title
  wrapped_title <- paste0(wrapped_title, title)

  return(wrapped_title)
}

Debugging

Individual Genes

• PCA with genes from signatures only

Violin Expression Plots

IndividualGenes_Violins(data = corrcounts_merge, metadata = metadata_merge, genes = c("CDKN1A", "CDKN2A", "GLB1","TP53","CCL2"), GroupingVariable = "Condition", plot=T, ncol=NULL, nrow=2, divide="CellType", invert_divide=FALSE,ColorValues=senescence_triggers_colors, pointSize=2, ColorVariable="SenescentType", title="Senescence", widthTitle=16,y_limits = NULL,legend_nrow=4, xlab="Condition",colorlab="") 

Using gene as id variables

Correlation Heatmap

CorrelationHeatmap(data=corrcounts_merge, 
                   metadata = metadata_merge, 
                   genes=c("CDKN1A", "CDKN2A", "GLB1","TP53","CCL2"), 
                   separate.by = "Condition", 
                   method = "pearson",  
                   colorlist = list(low = "#3F4193", mid = "#F9F4AE", high = "#B44141"),
                   limits_colorscale = c(-1,0,1), 
                   widthTitle = 16, 
                   title = "test", 
                   cluster_rows = TRUE, 
                   cluster_columns = TRUE,  
                   detailedresults = FALSE, 
                   legend_position="right",
                   titlesize=20)

Warning: Heatmap/annotation names are duplicated: pearson's coefficient
Warning: Heatmap/annotation names are duplicated: pearson's coefficient, pearson's coefficient
Warning: `legend_height` you specified is too small, use the default minimal height.
Warning: `legend_height` you specified is too small, use the default minimal height.
Warning: `legend_height` you specified is too small, use the default minimal height.

Expression Heatmaps

ROC/AUC

Cohen’s d

PCA with genes from signature only

#' @importFrom edgeR DGEList
#' @importFrom stats prcomp
#' @import ggplot2
#' @importFrom ggpubr ggarrange

plotPCA <- function(data, metadata, genes=NULL, scale=FALSE, center=TRUE, PCs=list(c(1,2)), ColorVariable=NULL,ColorValues=NULL,pointSize=5,legend_nrow=2, legend_position=c("bottom","top","right","left"),ncol=NULL, nrow=NULL){
  
  legend_position <- match.arg(legend_position)
  
  if (is.null(genes)){
    genes <-  row.names(data)
  }
  
  data <- data[row.names(data) %in% genes, , drop=F]
  
  if (!nrow(data)>1) stop(paste0("Error: Number of genes should be >1; In your data you have only found the gene ",genes))
  
  # Ensure metadata matches sample order if provided
  if (!is.null(metadata)) {
    colnames(metadata)[1] <- "Sample"
    rownames(metadata) <- metadata$Sample
    metadata <- metadata[colnames(data), , drop = FALSE]
    y <- edgeR::DGEList( log2(data+1), samples= metadata)  
  } else {
    y <- edgeR::DGEList( log2(data+1))  
  }
  

  
  nPCs <- max(unlist(PCs)) # get the maximum number of PC based on the user's choice
    
  PCAdata <- stats::prcomp(t(y$counts), scale=scale, center=center)
  PCAcounts <- PCAdata$x
  PCAcounts <- as.data.frame(PCAcounts)
  
  if (nPCs > ncol(PCAcounts)) stop("Error: Number of genes too low for number of chosen PCs. Please reduce number of PCs.")
      
  PCAcounts <-  cbind(PCAcounts[,1:nPCs],y$samples) 
  
  pltList <- list()
  
  for (pc in PCs){
    pc <- unlist(pc)
    ev = PCAdata$sdev^2
    pc_x <- round(100*ev[pc[1]]/sum(ev),2) 
    pc_y <- round(100*ev[pc[2]]/sum(ev),2) 
    
    plt <- ggplot2::ggplot(PCAcounts, ggplot2::aes_string(y = paste0("PC",pc[1]), x =  paste0("PC",pc[2])))
    
    
    # Add jittered points, optionally colored by ColorVariable.Default: Brewer Pallette "Paired"
    if (!is.null(ColorVariable)) {
      plt <- plt + ggplot2::geom_point(ggplot2::aes_string(fill = ColorVariable), size = pointSize, alpha = 0.5, shape=21, color="black")
    } else {
      plt <- plt + ggplot2::geom_point(size = pointSize, alpha = 0.5, shape=21, color="black", fill="#D8D8D8")  
    }
    
    # If ColorValues is provided, use a manual color scale; otherwise, if ColorVariable is provided,
    # use a default brewer palette.
    if (!is.null(ColorValues)) {
      plt <- plt + ggplot2::scale_fill_manual(values = ColorValues)
    } else if (!is.null(ColorVariable)) {
      plt <- plt + ggplot2::scale_fill_brewer(palette = "Paired")
    } 
    
    # Add axis labels (including variance)
    
    xlab <- paste0("PC",pc[1],": ",pc_x,"% variance")
    ylab <- paste0("PC",pc[2],": ",pc_y,"% variance")
    titleplot <- paste0("PC",pc[1], " vs PC",pc[2])
    
    plt <- plt + ggplot2::labs(fill = "", x = xlab, y = ylab, title=titleplot)
    
    # Change theme
    plt <- plt +
      ggplot2::theme_bw()+
    ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust=1),
          plot.title = ggplot2::element_text(hjust = 0.5),
          legend.position = "bottom")
    
    # Adjust legend rows if legend_nrow is specified
    if (!is.null(legend_nrow)) {
      plt <- plt + ggplot2::guides(fill = ggplot2::guide_legend(nrow = legend_nrow, position = legend_position))
    }
    
    # Add reference lines
    plt <- plt +
      ggplot2::geom_vline(xintercept=0, linetype="dotted") + 
      ggplot2::geom_hline(yintercept=0, linetype="dotted")  
    
    
    pltList <- c(pltList, list(plt))
    
    
    
  }
  
  n <- length(pltList)
  
  if(n==1){
    plt <- pltList[[1]]
  } else {
    
    # Determine grid layout
    if (is.null(ncol) && is.null(nrow)) {
      
      ncol <- ceiling(sqrt(n))
      nrow <- ceiling(n / ncol)
      
    } else if (is.null(ncol)){
      
      ncol <- ceiling(n / nrow)
      
    } else if (is.null(nrow)){
      
      nrow <- ceiling(n / ncol)
      
    }
    
    if (!is.null(ColorVariable)) {
      plt <- ggpubr::ggarrange(plotlist = pltList, ncol = ncol, nrow = nrow, common.legend = TRUE, align = "h")
    } else {
      plt <- ggpubr::ggarrange(plotlist = pltList, ncol = ncol, nrow = nrow, align = "h")
    }
    
    
    
  }
  
  print(plt)
  
  invisible(list(plt = plt,
              data = PCAcounts))
  
}